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We give a method to compute the smooth part of the density of states in a semi-classical expansion, 
when the Hamiltonian contains a Coulomb potential and non-cartesian coordinates are appropriate. 
We apply this method to the case of the hydrogen atom in a magnetic field with fixed z-component of 
the angular momentum. This is then compared with numerical results obtained by a high precision 
finite element approach. The agreement is excellent especially in the chaotic region of the spectrum. 
The need to go beyond the Thomas-Fermi model is clearly established. 



Qh' Studies in quantum chaos require a knowledge of the density of states of a given system in the semi-classical limit. 
This quantity is decomposed into two different parts, a smooth one and an oscillating one, related to periodic orbits 
■ of the corresponding classical system. A good knowledge of the smooth part is needed to unfold the spectrum before 
' determining its statistical behavior. Simple examples show how a bad unfolding deteriorates the statistics 

One of the standard tools to determine the smooth part of the density of states consists in computing its Laplace 
I""!' transform, i.e. the partition function in the semi-classical limit where a dimensionless h tends to zero but the tem- 
, perature remains fixed. Quite often the partition function is expressed by a functional integral of the Feynman-Kac 
' type 0]. But such an approach fails if the Hamiltonian contains a Coulomb potential. One way out of this difhculty 
is not to use cartesian coordinates, but more appropriate ones. Unfortunately it is difficult to work with a hmctional 
integral in non-cartesian coordinates despite progress made in some specific cases [H. 

We have therefore considered another approach based on the decomposition of the Hamiltonian into a differential op- 
erator and a multiplication operator. The differential operator is proportional to a dimensionless ft^ and its propagator 
is supposed to be known exactly in some non-cartesian coordinates. This formulation of the problem is particularly 
well adapted to the case where the multiplication operator contains a Coulomb potential, since in this case the use of 
semi-parabolic coordinates suppresses the Coulomb singularity at the origin. There is however a price to be paid in 
, the sense that the energy appears in the Hamiltonian. Nevertheless we have been able to compute the smooth part 
' of the density of states up to order . We have applied these formulas to the case of the hydrogen atom in a uniform 
magnetic field, which is a paradigmatic model in the study of quantum chaos 0, In this case we work in the fixed 
Lz vector. To the best of our knowledge these results arc new. 

We then compared the results with numerical energy levels obtained by means of a high precision finite element 
approach to the problem of the hydrogen atom in a uniform magnetic field. Details of the numerical method will 
be given elsewhere. The comparison of the analytical results with the numerical ones shows a very good agreement 
except in the lower part of the spectrum where the quasi-degeneracies specific to the hydrogen atom would require a 
Ch ' different treatment. The comparison demonstrates the need to go beyond the Thomas- Fermi type part of the density 
\ of states in order to satisfactorily account for the numerical results. 
• '~j ■ Our approach suggests that it should also be possible to derive the oscillating part of the density of states in a satis- 
/\ ' factory way. It is probably related to the periodic orbits in the way given in the literature 0] but to the best of our 
, knowledge, the usual derivation are not appropriate to the case of the hydrogen atom in a magnetic field. 
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I. SEMI-CLASSICAL EXPANSION OF THE SMOOTH PART OF THE DENSITY OF STATES: A 

GENERAL METHOD 



A problem of great interest, and which has often been treated consists in finding an asymptotic expansion for the 
density of states of an hamiltonian in terms of dimensionless parameter S corresponding to Planck constant. 

One way to achieve this is to find an asymptotic expansion for its Laplace transform, which is the partition function 
Z{t) given by 



Z{t) = Tr e' 



-tH 



(1) 



and to go back to the density of states by the inverse Laplace transform of Z{t). The following difficulty appears 
however in this program: t can be very large compared to for example. If we ignore this difficulty however by 
considering that t is fixed and 5 is small, we will get an expansion for the so-called smooth part of the density of 
states and will miss the oscillating part discovered by Balian-Bloch 0] and Gutzwiller Q ■ They are associated to 
terms like exp— in the expansion of Zit). This is the strategy we will follow. 
We will here consider the case where the hamiltonian is of the form 



and the propagator 



H = S^Ho + V 



Ut {r\r') = L-'i^'Ho+v)-] (^1^,) 



(2) 



(3) 



being a differential operator and V a multiplication operator. Our idea is the following. Consider that in the 
propagator given by Eq. [31 r is fixed and consequently so is V{r), we can then rewrite the propagator in the form 



where 

Wt (r|r') = L-*(«'^"+s)l {r\r') 

B being a multiplication operator defined by 

{B^j) (/) = [V{r') - V{r)] V'(r') 
We can now use the iterative representation of Wt (j'\r'), namely 

Wt {r\r') - Kt (r|r') + 

n=l •^*>*i 

where 



(4) 
(5) 
(6) 



[Kt-t,BKt,^t,B---Kt^]{r\r') (7) 



>t2>--->t„ 



Kt (r|r') 



(r|r') 



(8) 



In the standard case Hq = —A, one can check that the series (O corresponds to an asymptotic expansion in 5, when 
we use cartesian coordinates. The interest of this approach lies however in the fact that we can treat systems where 
more general coordinate systems are useful and particularly when the potential V contains a Coulomb potential, for 
which none of the standard techniques could work, to the best of our knowledge. 

II. POTENTIAL WITH A CYLINDRICAL SYMMETRY AND PARABOLIC COORDINATES 



We will consider the case of a three dimensional system with a potential V invariant under the rotation around the 
z-axis. We can therefore stay in the subspacc corresponding to ~ m. In what follows, m will be fixed, even when 
h — >■ 0. It is appropriate to use the following parabolic coordinates: 



X = ^Jxy cos 
Y = y^aJy sin 



(9) 



in which a Coulomb potential /ir ^ becomes 



+y' 

The dimensionless hamiltonian we will consider are of the form, 

H = -S^A, 



where the potential V has been decomposed into a Coulomb part and a more regular one V. 
The eigenvalue equation iJ?/' = -E'V' '^iU be written as 



where 



and 



(10) 





ip = /lip 


(11) 


J\ — ' fix 




[iZ) 


d d 
n-x = — tX—. — 
dx dx 


9 

m e 


(13) 


V{x,v)^[V- 


-Eo] 2 . 


(14) 



Here e = Eq — E > 0, Eq being the possible ionization threshold, which may depend on S. We now consider that 
equation pip quantizes ^ into eigenvalues ^ij{E) when e > 0. If we define the partition function 



g = Tr = ^ e-*^^ 



(15) 



and note that ~ — {tpj, ^-j^V'i) < which implies 



and finally the partition function 



Q{^l-^i,) = Q{E-E,{^J)). 



Hence the density of states we are looking at 

AA(i?,M) = ^e(£;-i?,(M)) 



(16) 



(17) 



(18) 



where fi is given, can be obtained once we know the partition function Q. 
But 



dxdyUt ixy\xy) 



(19) 



so that our task consists now in computing perturbatively G using equations (fTTj) . p2|) and ([T4|l . For this purpose we 
need to compute 



Kt {xy\x'y') 



{xy\x'y') 



(20) 



This can be done by using the fact that the eigenfunctions of S'^h^ are ipn{^^ where 



^^{x)=x\'^\e-^ L\T\x)^ 



[n + m)\ 



(21) 



E^^\x) being a Laguerre polynomial and 7 = The final result can be expressed as 



Kt ixy\x'y' 



^g-4|m|e 



AiTtS^ sinh 6 ^xx'yy' V tanh (9) 



exp 



' xx' + ■\/ yy' 



29 



I,r 



I xx' I. 



t52sinh(6l) ; Vi'5^sinh(6') 



261 



VV' (22) 



where 9 = and we introduced for later purpose the function 

= V2^e^"/™(x) (23) 
lm(x) being the usual Bessel function. Its usefulness comes from its simple asymptotic expansion 



'2 I I O ' 



exp 



kt [xy\x'y') 



[-^ [(V^-V^)' + (Vy-V^)'] -f + 



(24) 



The semi-classical limit of the propagator Kt that we will denote by Kt is obtained by taking 5 — > 0, \m\5 — > 
and (x,x' ,y,y') fixed 



(25) 



III. THE PARTITION FUNCTION 

We will compute Ut and therefore Q up to order (5". We decompose 



using the simplifying notation r = {x,y), 



Ut = + ul + uf 



(r|r') = e-*^M j^t (r|r) 



(26) 



(27) 



C/i (r|r') = e-*^W / Kt-tAr\r') KtAr'\r) f (r') - T>(r) 

Jt>ti>o 



dr' , 



(28) 



(r|r') = e-*^W / (rjn) ift,„t, (rajr) V{n)-V{r) V{r2)-V{r) 

Jt>ti>t2>0 

Let us begin by computing Ut to order (5*^. For this purpose we replace Kt by Kf 



Kt-t, {r\r')Kt, {r\r') 



exp -U 



' xx' + ^/yy' 



(47r)25Va;ya;'y'(t - tx)tx 



exp 



232 



where 



2t 



6^ it~h)h 



Let us now make the change of variables 



X' = \fx + 8u 



(29) 



(30) 



(31) 



(32) 
(33) 



with S being small then 

y(r') - V{r) = 2S \v^y/^u + Vyy/yv 



26' 



2" 



Vx — + Vxxxu +Vy— + Vyyyv + 2Vxy ^/xijuv 



where we have denoted by Vx = dxV. Introducing in addition 

x + y 



V = V + €- 



(34) 



it yields 



-tv 



dti 



271^5^^ Jo (t-ii)ti 



dudv exp 



VxVxu + Vyy/yv 



Vx-^ + VxxXu' + Y + ^yyV^'^ + '^Vxy^/xyuv 



dti 



so that finally one finds for Uf 



et 
"2 



VxX + Vyy]+-[Vx + Vy]+ VxxX + Vyyy 



Ulix.y) 



■J 



Vx (1 - etx) + Vy{l~ ety) + 2VxxX + 2Vyyy + 0{5) 



(35) 



(36) 



It is worth noting that we have replaced the domain of integration in the variable u: u > — by u > —00. This 
accounts to neglect terms of the order exp (— |). Such terms ( essential singularities in 6 ) arc responsible for the 
oscillating terms in the density of states, related to the classical periodic orbits. 

The second term U^{x, y) is treated in a similar way, replacing again Kt by Kt in equation (|29p . one finds 



u^{x,y) 



2 „-ty 



Vlx + K^y 



^^-K^xy 

These two terms give a contribution to the partition function, that we denote by 

G'^JJ dxdy [Ul{x,y) + U^{x,y)\ . 
It remains to compute the part of the partition function that we will denote by Q'^ 

G"= [[ dxdye-'^^^^y^Kt{xy\xy). 



Using eguation ip^ . it is given more explicitly by 



g^=C [['-^e.J-tVix,y) 



J J 



xy 



x + y 



— ^ (cosh (0) 1) j ^ j /|,n| [—^) 



(37) 



(38) 



(39) 



(40) 



where 



C = 



g-4e|m| ^ 

r-TTTT and 7 ~ — =. 

87r7sinh(6') ' V2e 



Using the identity 

I\rn\{a)i\„^\ib) = 1 + (/|,„|(a) - 1 + i\rn\ib) - l) + (/|rn|(a) - l) [l\rn\ib) - l) , 

we decompose into three different parts that we treat differently. A typical term will be of the form 



A= I ^ 



'|m| 



7sinh [6) 



- 1 



exp 



7 sinh (6) 



(cosh(0)-l) \g{x) 



(41) 



(42) 



We decompose A as follow: A ^ Ai + A2 



, , f°° dx f X \ 



cxp 



7 sinh (9) 



(cosh (6*) - 1) 



(43) 



But 



. , dx 

A-2^ 



'|m| 



7 sinh {6) 



exp 



7 sinh {9) 



(cosh(0)-l) (g(x)-g(O)). 



(44) 



Ai = g(0)v/27r7sinh {0) 



1 



sinh {9) 2 sinh 1 



(45) 



and 



where 



A. 



lim , ^ 

(5-fO 7 sinh [9) 



1 

~ " 4 



(46) 



e '^^^ [3(2;) - ,9( 



(47) 



In order to obtain this result, we have replaced /|m|(x) by its asymptotic behavior. In this way we obtain if 



1 \m\ r— et { y 1 
47rt(52 2tt6 tt V 48 



f f dxdy^_,y^^^, ^ 1 1 dx^^_,^,^y, ^^^^^^ ^^g^ 



JJ 



xy 



4-487r 



xy 



T_ 

32 



5(0) 



Stt 



et 



dxdy 



-tV(xy) 



W 2 1 
— TO 7 

Stt V 4 



dxdy 



-tV{xy) 



xy 



Vx+Vy 



(49) 



^0 ^ I^Le-^v-Co.o) 



(50) 



where g{x) =J^ f,-tV(x,y)-^y _ 

In deriving these expressions we have implicitly used some assumptions about the potential V{x, y). All the expressions 
given are correct provided 



and a + i > 0. 



V{x, y) - V{0, y) - when x — > 
V{x, y) - V{x, 0) ~ 2/" when y — ^ 



IV. THE INTEGRATED DENSITY OF STATES 



(51) 
(52) 



The integrated density of states (IDOS) M{E) = J^j 'S) {E — Ej) can now be obtained from the knowledge of the 
partition function. But we must not forget that the ionization threshold E^ may also depend on 5. We assume that 
i?o = 5ei and consider that semi-classically the IDOS should be computed by keeping e = E — Eq fixed and taking S 
small. 



We therefore take 



V{x,y) = W{x,y)^5e^^^ (53) 



where W{x, y) = y) + e). We decompose now the IDOS into three parts corresponding to their importance 



2 

in the semi-classical limit 



Using the representation 



Af=^ + f + N,. (54) 



le-'^ = Jdt,e-''^^[^,-a]^ (55) 



we see that 



If we go back to the original coordinates (p, z) instead of the parabolic ones, we can rewrite this term in a more 
physical expression 

^fo = (^) ' / '^'P J dpdz e (^-e - ^ - V(p, r) + ^) (57) 

which is the Thomas-Fermi form of the density of states, except that the volume element is dpdz, not pdpdz and — e 
replaces E. In order to compute the other terms, we can use the identities 



We then found that 



M = -^V2^//^e(,-iy (.,,)) 

27r J J y/xy 



since W{0, y) = W{x, 0) = by our assumptions and 



ff'^'yei,~Wix,y))ix + y)-\^ (60) 



M, = ^[9m'^l^+^JJ^6{fi- Wix, y)) g{x, y) (61) 



where 



g{x,y) = em - jmlj -er 

' ^ ^ ' \vi^. + W„) 



,22 



TO — 



^ {xW^. + yW^yy) (62) 




FIG. 1: Plots of numerical energy counting versus the semi-classical integrated density of states for 5 = 4.64 ■ 10^^ (7 — 10"*), 
once the Thomas- Fermi approximation Mo{^i) and then with corrections up to second order M{ei). 



V. THE INTEGRATED DENSITY OF STATES OF THE HYDROGEN ATOM IN A UNIFORM 

MAGNETIC FIELD 

The case of the hydrogen atom in a magnetic field constant and directed along the z axis represents the most 
important application of our general formula. 

In term of the units of length aq = -^-^ and L = ^ -g^, our parameter 5 will be (5 = yj^- Here M is the 

2 

reduced mass, B the magnetic field strength and e the charge of the electron. In miits of energy the dimensionless 
hamiltonian reads 

A 1 9 1 



the parameter e is given by 



A777 

e^ — +5ei~E (64) 



where ei = ^(1 + \'m\) and finally the diamagnetic potential in this case is expressed as 

V(x,2/) = ^. (65) 

We may note that many authors 0, i, H, i, Ell, E [H, Q [H, [H, [13, El, [H, H IMi have introduced instead of 
our parameter 5 a parameter 7 = 5^ = ^clq. 

Let us use the parametrization c S [1, 00] instead of e S [0, 00] given by the relation 

2e=(c-l)c-i, (66) 
in this case the various terms in the IDOS become 



c '^T^ Jo Vy \ c-l + yj 
7Vi(e) = -|m|^= / ^ 1 + -— ^==: -(1 + |m|)c ^ / — ^ 

■^2(^) = K9-'-lVr^^' (67) 
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FIG. 2: Difference between the semi-classical integrated density of states for S — 4.64 • 10^^ (7 — 10~*) and its best fit curve 
y — + h computed in the non-shadowed region i £ [90; 300]. h is equal to 8.56 in the case of first order corrections only 
and 8.39 in the case of both first and second order corrections. The existence of a quasi-invariant at low energies is clearly 
distinguished by the vertical clustering produced by similar energy values. 



where we have defined for ease the functions 



(68) 



and 



liv) = a 



•1 + 



c-l + y 



+ asy + 



(C - l + y)2 



(69) 



together with the coefficients 




(70) 



(71) 



_1_ fm^ _l\ 
2^ \~Y 3j ' 



(72) 




(73) 



All these integrals are of the elliptic type and thus evaluations have been performed numerically. 
It is important to notice that the Thomas- Fermi IDOSs remains finite at the ionization threshold (e = 0). This is 
however not compatible with the diverging IDOS at this threshold as it corresponds to an accumulation point for the 
eigencnergies. We therefore expect that the corrections A/i(e) and A/'2(e) become more and more relevant in the limit 
e ^ 0. This will be confirmed in the next section by the numerical evaluation of the IDOS. 



The preceding formulas have been compared with numerical results of the energy levels in the — 0+ subspace of 
the hydrogen atom in a magnetic field strength of 23.5 T that is (5 = 4.64 • 10"^ (7 = 10~^). The latter were obtained 
through a high precision finite element approach in cylindrical coordinates. Details of the numerical method and its 
comparison with other approaches and results will be given in forthcoming paper [2^ . 

Choice has been made to present the results as a parametric function {JV{ei);i) where correspond to the i*'' 
numerical energies. Nevertheless with regards to the subspace considered, only even eigenfunctions are taken into 
account. Therefore we expect a linear correspondence with a 1/2 slope since, on average between each even energy 
lies an odd one. 

Figure [1] illustrates the parametric functions obtained through the Thomas-Fermi IDOS and its first and second order 
corrections with regards to S. Both functions exhibit a clear linear growth as expected. However there are increasing 
deviations of the Thomas- Fermi curve at high energies corroborating our last remarks in section |Vl namely that the 
Thomas-Fermi term is not valid as one approaches the ionization threshold. 

A magnification of the linear correspondance to the theoretical half slope is presented in Figure [2] where deviations 
from the best fit line y = ^x + h are reported once for the function A/'(e) and once without the second order correction 
i.e. M{e) — A2(£). h is left as a free parameter since its value cannot be inferred rigorously in the analytical approach 
(no absolute counting from below). As expected, one notices failure to match a linear approximation in the shaded 
regions corresponding to the low and very high energy range. Strong deviations are indeed observed in the former 
domain due to the existence of a quasi-invariant which generates quasi-degeneracies. 

The magnetic field is here not strong enough to take over the Coulomb term and one observes clearly on Figure [2] a 
clustering of the eigenvalues around the quasi-degeneracies for which, at the lower end, simple perturbation theory in 
the magnetic field strength would already provide an accurate picture. At the higher end of the lowest shaded region 
one would probably benefit from an analytical expression of the hamiltonian in the Gay-Delande's basis [lo| . and then 
apply a semi-classical expansion in order to account for the density of states. 

The deviations at high energies can be attributed to two presumed causes, namely missing higher correction orders 
or the departure from the semi-classical domain. However further analysis would be required in order to determine 
which of them is the first dominant. 

Although in the mid domain deviations from the theoretical fit line are small, it is not immediately clear how one 
could define the best domain boundaries for each curve. Moreover those boundaries are required in order to compute 
the h parameter. In order to resolve such ambiquities we plot in Figure [3] the RMS of residuals around linearity 
y = + h with h left as a free parameter for every possible intervals e^, k G Contributions arising from the 



VI. COMPARISON WITH NUMERICAL RESULTS 
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i 

FIG. 3: Contourplot for S = 4.64 ■ 10"^ (7 = 10"") of the RMS of residuals between the semi-classical integrated density of 
states and the theoretical curve given hy y = + h where the height h is left as a free parameter within the ek interval given 
by fc G The black curve corresponds to a calculation with first order correction only, the gray one with both first and 

second order corrections. 

first and second order corrections have been computed separately so that one can easily pictures the difference. The 
results show indeed that including first and second order corrections improves the overall linearity in different ways. 
One can think first at fixed interval where the RMS of residual would then be lower or second at constant RMS 
residual where the range would then be wider. The latter is particularly important as it clearly demonstrates the tiny 
refinement induced by the second order corrections at high energies. 

Considering the non-shadowed domain of Figure [21 the RMS of residuals equals 1.06 with all corrective terms and 
1.35 with just the first order correction. 

It is now quite interesting to compare the different contributions of Af{e) represented by Ao, Mi and A2 in Eq. ([M]) 
to the slope of figure [TJ For this purpose, the energy interval is taken from the preceding considerations, namely ek, 
k G [90,280]. The exact numerical values of the best fit line are presented in table [H It appears clearly again that 
the Thomas- Fermi term i.e. S~'^Afo, is definitively not able to reproduce quantitavely the counting function. A very 
important contribution is indeed brought by the first order correction Afi in terms of accuracy with respect to the 
theoretical slope value. At last, the second order correction Af2 further improves the slope value but with higher RMS 
of residuals. The lower RMS of residuals value when ignoring A/2 is an artefact since it was measured with respect to 
a less correct slope. Figure [31 gives a more consistant picture on the issue of analyzing the RMS fluctuations. In fact, 
we would like to stress that fluctuations of the energy level spacings are an important feature of the spectrum of the 
hydro gen atom under constant magnetic field, and carry out important statistical informations related to quantum 
chaos [iTI. [21I. [23l. [2^. [25l. [26l. [27| .Thcv will be analyzed in a forthcoming paper (2^ since these are out of the scope of 
the present paper. 

Finally a last remark: it should be clear to the reader that both first and second order corrections do depend on the 
magnetic quantum number m. Although we have only computed the case m = 0, corrections to the Thomas-Fermi 
term would be higher for m ^ according to Eqs. (|67ll70p . 



function 



slope height rms of residuals 




0.98 



0.48 



0.54 



TABLE I: Best linear least square approximations of the smoothed density of states in the energy range starting at the 90*'' 
even energy levels up to the 280*''. 
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